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Abstract 

This project work report provides a full solution of 
simplified Navier Stokes equations for The Incom- 
pressible Couette Problem. The well known analyt- 
ical solution to the problem of incompressible cou- 
ette is compared with a numerical solution. In that 
paper, I will provide a full solution with simple C 
code instead of MatLab or Fortran codes, which are 
known. For discrete problem formulation, implicit 
Crank-Nicolson method was used. Finally, the sys- 
tem of equation (tridiagonal) is solved with both 
Thomas and simple Gauss Method. Results of both 
methods are compared. 



1 Introduction 

Main problem is shown in figure (Q. There is 
viscous flow between two parallel plates. Upper 
plate is moving in x direction with constans velocity 
(U = Ue). Lower one is not moving ([/ = 0). We 
are looking for a solution to describe velocity vector 
field in the model (between two plates). 



* Thanks for Grzegorz Juraszek (for English languague 
checking) . 
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Figure 1: Schematic representation of Couette Flow 
Problem 

2 Fundamental Equations 

Most of incompressible fluid mechanics (dynam- 
ics) problems are described by simple Navier-Stokes 
equation for incompressible fluid velocity, which can 
be written with a form: 

— = {-iIS7)u-Vip + v\7'^u + g, (1) 

where tp is defined is defined as the relation of 
pressure to density: 



P 

P 



(2) 



and u is a kinematics viscosity of the fluid. 
We will also use a continuity equation, which can 
be written as follows ^: 



assume constans density of the fluid. 



D = V-v = (3) 

Of course, in a case of couette incompressible flow 
we will use several simplifications of ((T|l. 



3 Mathematical Formulation 
of the Couette Problem 

Incompressible couette problem is not needed to 
solve full Navier-Stokes equations. There is no ex- 
ternal force, so first simplification of will be: 



9u , ^ 
- = (-.V).- 



(4) 



In m there can be found easy proof that in cou- 
ette problem there are no pressure gradients, which 
means that: 



(5) 



We will ignore a convection effects so, equation 
Q can be written with a form: 



— = uV^u 
dt 



(6) 



Now we have simple differential equation for ve- 
locity vector field. That equation is a vector type 
and can be simplified even more. Let us write 
continuity equation Q in differential form. Let 
u — {u,v), then continuity equation can be ex- 
panded as follows^: 



du dv ^ 
dx dy 



(7) 



We know that there is no u velocity component 
gradient along x axis (symmetry of the problem). 



du 
dx 



= 



(8) 



Evaluation of Taylor series'^ at points y = and 
y = D gives us a proof that only one possible and 
physically correct value for y component of velocity 
u is: 

v = (9) 

Because of © equation © can be written as fol- 
lows: 



^Only with assumption of non compressible fluid. 

^ Those simple expressions can be found in 1 , chapter 



9.2. 



du 
'dt 



d^U 

dy'^ 



(10) 



Our problem is now simplified to mathematical 
problem of solving equations like H10() . That is now 
a governing problem of incompressible couette flow 
analysis. 

3.1 Analytical soulution 

Analytical solution for velocity profile of steady 
flow, without time-changes (steady state) can be 
found very easily in an equation: 



V—-0 
dy"^ 



(11) 



And without any changes of viscosity v it can be 
written in form: 



= 



(12) 



After simple two times integration of equation 
(|12|l we have analytical solution function of H12I) : 

u = ci-y + C2 (13) 
where ci and C2 are integration constans. 

3.2 Boundary Conditions for The 
Analytical Solution 

Simple boundary conditions are provided in that 
problem. We know that: 
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y=0 
y-D 



(14) 



Simple applying it to our solution H13|l gives a 
more specified one, where ci — ^ and C2 = 0: 



Ue 

u = — - y 



D 



(15) 



It means, that a relationship between u and y is 
linear. A Better idea to write that with mathemat- 
ical expression is: 



D 



(16) 



Where ^ is a constans for the problem (initial u 
velocity vs size of the model). 
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4 Numerical Solution 



4.1 Non-dimensional Form 

Let us define some new non-dimensional variables'*: 



V 

y = D 
t = 



(17) 



Now let us place these variables into equation 
and we have now a non-dimensional equation 
written as follows: 



du/i 



Now we replace all the variables to nondimen- 
sional, like defined in p7|l : 



P 



du 

W \D 



92 



(19) 



Now we will remove all chars from that equa- 
tion (only for simplification of notation), and the 
equation becomes^ to: 



du 

'at 



u 



(20) 



\DpueJ dy"^ 
In that equation Reynold's number Re appears, 
and is defined as: 



Re = 



DpUe 



(21) 



Where Re is Reynold's number that depends on 
D height of couette model. Finally, the last form of 
the equation for the couette problem can be written 
as follows: 



(22) 



du _ f dhi 
dt ^ \Re) 9y2 

We will try to formulate numerical solution of the 
equation if^ . 



4.2 Finite - Difference Representa- 
tion 

In our solution of equation H22I) we will use Crank- 
Nicolson technique, so discrete representation of 
that equation can be written ^ as: 

^Exacly the same, like in T'. 
^Constans simplification also implemented 
^That representation is based on central discrete differen- 
tial operators. 



2{Ay)^Re' 



iu]tl+u]+^-2u]+^+ 



n-|-l I n-l-1 



(23) 

Simple grouping of all terms which are placed in 
time step (n-|-l) on the left side and rest of them 
- on right side, gives us an equation which can be 
written as: 



-1) 



Ai 



Bu 



n+l 



(24) 



Where Kj is known and depends only on values 
u at n time step: 



K - [I 
' ^ {AyfRe 



At 



^ 2{AyYRe^ ' 



(25) 



Constans A and B are defined as follows^: 
A 



B = 1 



2{AyYRe 

At 
{AyYRe 



(26) 
(27) 



4.3 Crank - Nicolson Implicit scheme 

For numerical solution we will use one-dimensional 
grid points (1, . . . , TV -|- 1) where we will keep cal- 
culated u velocities. That means u has values from 
the range (mi,M2, . . . ,un+i). We know (from fixed 
boundary conditions), that: mi = and mat 4.1 — 0. 
Simple analysis of the equation (|24|l gives us a sys- 
tem of equations, which can be described by matrix 
equation: 



(28) 



Where A is tridiagonal [A^ — 1] • [A^ — 1] matrix 
of constant A and B values: 



B A 
ABA 



ABA 
A B 



X vector is a vector of u values: 



(29) 



X = [ui,U2, ■ 



(30) 



^Directly from equation 1231 . 
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Y vector is a vector of constans Kj values: 



Y =[Ki,...,Kn+i] 



(31) 



5 Solving The System of Lin- 
ear Equations 

Now the problem comes to solving the matrix - vec- 
tor equation H28(l . There are a lot of numerical 
methods for that®, and we will try to choose two 
of them: Thomas and Gauss method. Both are 
very similar, and I will start with a description of 
my implementation with the simple Gauss method. 

5.1 Gauss Method 

Choice of the Gauss method for solving system of 
linear equations is the easiest way. This simple al- 
gorithm is well known, and wc can do it very easily 
by hand on the paper. However, for big matrices 
(big N value) a computer program will provide us 
with a fast and precise solution. A very important 
thing is that time spent on writing (or implement- 
ing, if Gauss procedure was written before) is very 
short, because of its simplicity. 

I used a Gauss procedure with partial choice of 
a/the general element. That is a well known tech- 
nique for taking the first element from a column of 
a matrix for better numerical accuracy. 

The whole Gauss procedure of solving a system of 
equations contains three steps. First, we are look- 
ing for the general element. 

After that, when a general element is in the first 
row (we make an exchange of rows^) we make some 
simple calculations (for every value in every row 
and column of the matrix) for the simplified matrix 
to be diagonal (instead of a tridiagonal one which 
we have at the beginning). That is all, because af- 
ter diagonalization I implement a simple procedure 
(from the end row to the start row of the matrix) 
which calculates the whole vector X. There are my 
values of Ui velocity in all the model^°. 

5.2 Thomas Method 

Thomas' method, described in pP is simplified ver- 
sion of Gauss method, created especially for tridiag- 



* Especially for tridiagonal matrices like A 
^We made it for A matrix, and for X, Y too. 
^"^More detailed description of Gauss method can be found 
in a lot of books on numerical methods, like 2]. 



onal matrices. There is one disadvantage of Gauss 
method which disappears when Thomas' method is 
implemented. Gauss method is rather slow, and 
lot of computational time is lost, because of special 
type of matrix. Tridiagonal matrices contain a lot 
of free (zero) values. In the Gauss method these 
values joins the calculation, what is useless. 

Thomas' simplification for tridiagonal matrices is 
to get only values from non-zero tridiagonal part of 
matrix. Simply applying a Thomas' equations for 
our governing matrix equation (|28(l gives us: 



d, = B 



Ui = Ui 



u,A 



d- 



(32) 



(33) 



We know that exact value of um is defined as 
follows: 



UM 



"-M 

B 



(34) 



Now solution of the system of equations will be 
rather easy. We will use recursion like that: 



B 



(35) 



That easy recursion provides us a solution for the 
linear system of equations. 

6 Results 

Main results are provided as plots of the function: 



D 



= .f 



(36) 



In figure (01 there is drawn an analytical solu- 
tion to the problem of couette flow. That is linear 
function, and wc expect that after a several time 
steps of numerical procedure we will have the same 
configuration of velocity field. 

6.1 Different Time Steps 

In figure ||2J) there are results of velocity u calcu- 
lation for several different time steps. Analytical 
solution is also drawn there. 

As we can see in the figure (jSJ - the solution is 
going to be same as analytical one. Beginning state 
(known from boundary conditions) is changing and 
relaxing. 
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Figure 2: Analytical exact solution. 
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Figure 3: Results for different time steps of numer- 
ical calculations. 

6.2 Results for Different Reynolds 
Numbers 

In the figure Q there is plot of numerical calcula- 
tions for different Reynold's numbers. For example 
Reynold's number depends on i.e. viscosity of the 
fluid, size of couette model. As it is shown on the 
plot there is strong relationship between the speed 
of the velocity field changes and Reynold's number. 
In a couple of words: when Reynolds number in- 
creases - frequency of changes also increases. 

6.3 Results for Different Grid Den- 
sity 

In figure lO there is an example of calculations of 
velocity field for different grid density (TV number). 



Figure 4: Calculation for same time (t — 2500) and 
different Reynold's numbers. 



We see that there is also strong correlation between 
grid density, and speed of changes on the grid. Also, 
very interesting case = 10 shows, that for low 
density of the grid changes are very fast, and not 
accurate. 




Figure 5: Calculation for same time {t = 2500), 
same Reynold's numbers (Re=5000) and different 
grid density (A^ number). 



6.4 Conclusion 

Solving of Incompressible Couette Problem can be 
good way to check numerical method, because of 
existing analytical solution. In that report there 
were presented two methods of solving system of 
equations: Gauss and Thomas' method. System of 
equations was taken from Crank-Nicolson Implicit 
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scheme. Well known linear relationships were 
served. 



7 APPENDIX A 



#include <stdlib.h> 
#include <stdio.h> 
#include <iiiath.h> 

#define N (40) 
#defiiie NN (N+1) 

void Zamien (double *a, double *b) { 
double c; 

c=*a; *a=*b; *b=c; 

} 

void WypiszMacierz (double A[NN] [NN] , int n) { 
int i , j ; 

for(j=l;j<n;j++) 
{ 

for(i=l;i<n;i++) 
i 

printf ("7,2.4f 

} 

printf ("\n") ; 

} 

} 

void Gauss (double A[NN] [NN] , double *b, double *x, int n) { 
int i, j ,k; 
double m; 

// Gauss Elimination 

for(i=0;i<n;i++) 

// Step #1 : Change governing element 

m=fabs(A[i] [i] ) ; 
k=i; 

for(j=i+l;j<n;j++) 
if (fabs(A[i] [j])>m) 
i 

m=fabs(A[i] [j]) ; 
k=j; 

} 

if (k!=i) 

for(j=0;j<n;j++) 
i 

Zamien (&A [ j ] [i] , &A [ j ] [k] ) ; 
Zamien (&b[i+l] ,&b[k+l]); 



// show matrix 
",A[i] [j]); 
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} 



// Step #2: make it tricmgle 

for(j=i+l;j<n;j++) 

{ 

m = A[i] [j]/A[i] [i] ; 

for(k=i;k<ii;k++) 

A[k] [j] = A[k] [j] - m*A[k] [i] ; 

b[j+l] = b[j+l] - m*b[i+l]; 

} 

} 

// Step#3: Solve now 

for(i=n-l;i>=l;i — ) 
{ 

for(j=i+l; j<ii; j++) 

b[i+l] = b[i+l]-A[j] [i]*x[j+l] ; 
x[i+l] = b[i+l]/A[i] [i] ; 

} 



int main (void) 
{ 

double U [N*2+2] ={0} , A [N*2+2] ={0} , B [M*2+2] ={0} , C [M*2+2] ={0} , D [N*2+2] ={0} , Y [M*2+2] ={0} ; 
// initialization 

double OneOverN = 1 .0/ (double) N; 

// Reynolds number 
// dt parameter 



double Re=5000; 
double EE=1.0; 
double t=0; 

double dt=EE*Re*0ne0verN*2; 

double AA=-0.5*EE; 
double BB=1.0+EE; 
int KKEND=1122; 
int KKM0D=1; 
int KK; 
int i, j ,k; 
int M; 

double GMatrix [NN] [NN] ={0} ; 
double test; 
Y[1]=0; // init 



// delta time 



// for a loop 
// for loops too 
// temporary needed variable 
// for Gauss Elimination 



// apply boundary conditions for Couette Problem 

U[1]=0.0; 
U[NN]=1.0; 



// initial conditions (zero as values of vertical velocity inside of the couette model) 

for(j=2;j<=N;j++) 
U[j]=0.0; 
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A[1]=B[1]=C[1]=D[1]=1.0; 

f or ( KK= 1 ; KK<=KKEND ; KK++ ) 

for(j=2;j<=N;j++) 
i 

Y[j]=Y[j-l]+OneOverN; 

A[j]=AA; 

if (j==N) 

A[j]=0.0; 
D[j]=BB; 
B[j]=AA; 
if(j==2) 

B[j]=0.0; 

C [ j ] = ( 1 . 0-EE) *U [ j ] +0 . 5*EE* (U [ j +1] +U [j - 1] ) ; 
if (j==N) 

C[j]=C[j]-AA*U[NN] ; 



// Gauss 
// C[] 

// A[]B[]D[] 
// U[] 



free 

for matrix calculation 
X 



// calculate matrix for Gauss Elimination 

GMatrixCO] [0]=D[1] ; 
GMatrix[l] [0]=A[1] ; 



for(i=l;i<N-l;i++) 
i 

GMatrix[i-l] [i]=B[i+l] 
GMatrixCi] [i]=D[i+l] ; 
GMatrix[i+l] [i]=A[i+l] 

> 



// GMatrix[l] [2] =B [2] 
// GMatrix[2] [2]=D[2] 
// GMatrix [3] [2] =A [2] 



GMatrix[N-2] [N-1]=B[N] 
GMatrix [N-1] [N-1]=D[M] 



Gauss (GMatrix, C,U,N) : 



// Gauss solving function 



Y[1]=0.0; 

Y [NN] =Y [N] +OneOverN ; 



t=t+dt ; 

test=KK '/, KKMOD; 



// time increment 



if (test < 0.01) 
i 

printf ("KK,TIME\n") ; 
printf ("7.d,7.f\n",KK,t) ; 



// print the results 
// info 1 
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printf (",J,Y[J] ,U[j]\n"); // info 2 
for(j=l; j<=NN; j++) 

printf Cy.d , y.f, y.f\n",j,U[j] ,Y[j]); 

printf ("\n \n \n \n"); // for nice view of several datas 

} 

} 

return (1) ; 

> 
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#include <stdlib.h> 
#include <stdio.h> 

#define N (50) 
#define NN (N+1) 

int main (void) 
{ 

double U[N*2+2] ,A[N*2+2] ,B[N*2+2] ,C[N*2+2] ,D[N*2+2] ,Y[N*2+2] ; 
// initialization 

double OneOverN = 1.0/ (double) N; 
double Re=7000; 
double EE=1.0; 
double t=0; 

double dt=EE*Re*0ne0verN*2; 
double AA=-0.5*EE; 
double BB=1.0+EE; 
int KKEND=122; 
int KKM0D=1; 
int KK; 
int j ,k; 
int M; 

double test ; 
Y[1]=0; // init 

// apply boundary conditions for Couette Problem 
U[1]=0.0; 
U[NN]=1.0; 

// initial conditions (zero as values of vertical velocity inside of the couette model) 
for(j=2;j<=N;j++) 

U[j]=0.0; 
A[1]=B[1]=C[1]=D[1]=1.0; 

printf ("dt=7.f, Re=y.f , N=7.d \n",dt,Re, N) ; 



// Reynolds number 
// dt parameter 

// delta time 



// for a loop 

// for loops too 

// temporary needed variable 
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for (KK=1 ; KK<=KKEND ; KK++) 
{ 

for(j=2;j<=N;j++) 
i 

Y[j]=Y[j-l]+OneOverN; 

A[j]=AA; 

if (j==N) 

A[j]=0.0; 

D[j]=BB; 
B[j]=AA; 

if(j==2) 

B[j]=0.0; 

C [ j ] = ( 1 . 0-EE) *U [ j ] +0 . 5*EE* (U [ j +1] +U [j - 1] ) ; 

if (j==N) 

C[j]=C[j]-AA*U[NN] ; 



// upper bidiagonal form 

for(j=3;j<=N;j++) 

{ 

D[j]=D[j]-B[j]*A[j-l]/D[j-l]; 
C[j]=C[j]-C[j-l]*B[j]/D[j-l]; 

> 



// calculation of U[j] 

for(k=2;k<N;k++) 

i 

M=N-(k-2) ; 

U [M] = (C [M] -A [M] *U [M+1] ) /D [M] 

} 



// Appendix A 



Y[1]=0.0; 

Y [NN] =Y [N] +OneOverN ; 
t=t+dt ; 

test=KK '/, KKMOD; 



// time increment 



ifCtest < 0.01) 
i 

printf ("KK,TIME\n") ; 
printf ("•/,d//,f\n",KK,t); 



// print the results 
// info 1 



printf (",J,Y[J] ,U[j]\n"); // info 2 

for(j=l;j<=NN;j++) 

printf ("7.d , y.f, y.f\n",j,U[j] ,Y[j]); 

printf ("\n \n \n \n"); // for nice view of several datas 



return (1) ; 
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